#' ---
#' title: "Polar Maps in Leaflet"
#' author: "Bhaskar V. Karambelkar"
#' ---

library(leaflet)

#' ## Artic Projections

#' There is a [polarmap.js](http://webmap.arcticconnect.org/)
#' leaflet plugin available, but that one is not easy to integrate in to the R package.<br/>
#' But thankfully it does provide Tiles in different projections
#' which can be used with Proj4Leaflet.
#' In all it supports 6 projections and corresponding tile layers.
#' <br/>
#' The polarmap.js supports only Artic data, for Antartica see the end of this document.


#' All these numbers and calculations come from the polarmap.js plugin, specifically from these files
#'
#' - http://webmap.arcticconnect.org/polarmap.js/dist/polarmap-src.js
#' - http://webmap.arcticconnect.org/tiles.html
#' - http://webmap.arcticconnect.org/usage.html
#'
extent <- 11000000 + 9036842.762 + 667
origin = c(-extent, extent)
maxResolution <- ((extent - -extent) / 256)
defZoom <- 4
bounds <- list(c(-extent, extent),c(extent, -extent))
minZoom <- 0
maxZoom <- 18
resolutions <- purrr::map_dbl(minZoom:maxZoom,function(x) maxResolution/(2^x))

# 6 Projection EPSG Codes
projections <- c('3571', '3572', '3573', '3574', '3575', '3576')
# Corresponding proj4defs codes for each projection
proj4defs <- list(
  '3571' = '+proj=laea +lat_0=90 +lon_0=180 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs',
  '3572' = '+proj=laea +lat_0=90 +lon_0=-150 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs',
  '3573' = '+proj=laea +lat_0=90 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs',
  '3574' = '+proj=laea +lat_0=90 +lon_0=-40 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs',
  '3575' = '+proj=laea +lat_0=90 +lon_0=10 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs',
  '3576' = '+proj=laea +lat_0=90 +lon_0=90 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs'
)

# create a CRS instance for each projection
crses <- purrr::map(projections, function(code) {
  leafletCRS(
    crsClass = 'L.Proj.CRS',
    code = sprintf("EPSG:%s",code),
    proj4def = proj4defs[[code]],
    origin = origin,
    resolutions = resolutions,
    bounds = bounds
  )
})

# Tile URL Template for each projection
tileURLtemplates <- purrr::map(projections, function(code) {
  sprintf('http://{s}.tiles.arcticconnect.org/osm_%s/{z}/{x}/{y}.png',
          code)
})

# We can't add all 6 tiles to our leaflet map,
# because each one is in a different projection,
# and you can have only one projection per map in Leaflet.
# So we create 6 maps.
polarmaps <- purrr::map2(crses, tileURLtemplates,
    function(crs, tileURLTemplate) {
      leaflet(options= leafletOptions(
        crs=crs, minZoom = minZoom, maxZoom = maxZoom)) %>%
        setView(0, 90, defZoom) %>%
        addTiles(urlTemplate = tileURLTemplate,
          attribution = "Map © ArcticConnect. Data © OpenStreetMap contributors",
          options = tileOptions(subdomains = "abc", noWrap = TRUE,
                      continuousWorld = FALSE))
    })

#' #### EPSG:3571
polarmaps[[1]] %>%
  addGraticule()

#' #### EPSG:3572
polarmaps[[2]]

#' #### EPSG:3573
polarmaps[[3]]

#' #### EPSG:3574
polarmaps[[4]]

#' #### EPSG:3575
polarmaps[[5]]

#' #### EPSG:3576
polarmaps[[6]]

#' ## Antartica
#' Code adapted from
#' https://github.com/nasa-gibs/gibs-web-examples/blob/release/examples/leaflet/antarctic-epsg3031.js <br/>

resolutions <- c(8192, 4096, 2048, 1024, 512, 256)
zoom <- 0
maxZoom <- 5

border <- geojsonio::geojson_read(system.file('examples/Seamask_medium_res_polygon.kml', package = 'leaflet'), what='sp')
points <-  geojsonio::geojson_read(system.file('examples/Historic_sites_and_monuments_point.kml', package='leaflet'), what='sp')

crsAntartica <-  leafletCRS(
  crsClass = 'L.Proj.CRS',
  code = 'EPSG:3031',
  proj4def = '+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs',
  resolutions = resolutions,
  origin = c(-4194304, 4194304),
  bounds =  list( c(-4194304, -4194304), c(4194304, 4194304) )
)

antarticaTilesURL <- "//map1{s}.vis.earthdata.nasa.gov/wmts-antarctic/MODIS_Aqua_CorrectedReflectance_TrueColor/default/2014-12-01/EPSG3031_250m/{z}/{y}/{x}.jpg"

leaflet(options= leafletOptions(
  crs=crsAntartica, minZoom = zoom, maxZoom=maxZoom, worldCopyJump = FALSE)) %>%
  setView(0, -90, 0) %>%
  addPolygons(data=border, color = '#ff0000', weight = 2, fill = FALSE) %>%
  addCircleMarkers(data=points, label=~Name) %>%
  addTiles(urlTemplate = antarticaTilesURL,
           layerId = "antartica_tiles",
           attribution = "<a href='https://earthdata.nasa.gov/gibs'> NASA EOSDIS GIBS</a>&nbsp;&nbsp;&nbsp; <a href='https://github.com/nasa-gibs/web-examples/blob/release/leaflet/js/antarctic-epsg3031.js'> View Source </a>",
           options = tileOptions(
             tileSize =512,
             subdomains = "abc",
             noWrap = TRUE,
             continuousWorld = TRUE,
             format = "image%2Fjpeg"
           )) %>%
  addGraticule() %>%
  htmlwidgets::onRender(
    "function(el, t){
       var myMap = this;
       debugger;
       var tileLayer = myMap.layerManager._byLayerId['tile\\nantartica_tiles'];

       // HACK: BEGIN
       // Leaflet does not yet handle these kind of projections nicely. Monkey
       // patch the getTileUrl function to ensure requests are within
       // tile matrix set boundaries.
       var superGetTileUrl = tileLayer.getTileUrl;

       tileLayer.getTileUrl = function(coords) {
         debugger;
         var max = Math.pow(2, tileLayer._getZoomForUrl() + 1);
         if ( coords.x < 0 ) { return ''; }
         if ( coords.y < 0 ) { return ''; }
         if ( coords.x >= max ) { return ''; }
         if ( coords.y >= max ) { return ''; }
           return superGetTileUrl.call(tileLayer, coords);
       };
       // HACK: END
    }")
